%%%%%%%%%%%%%%%%%%%%%%%%%
% explicit Euler's method
%%%%%%%%%%%%%%%%%%%%%%%%%
%
% consider the equation y'=y+2x, y(0)=-1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all
close all

  % A) input data h, N, x0, y0, f(x,y)
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

h = 0.1;           % step-size
N = 10;            % number of points
x0 = 0; y0 = -1;   % initial condition - starting point
plot(x0, y0, 'r*') % plot the starting point as a red asterisk
hold on            % continue on the same picture

  % explicit Euler
x = x0; y = y0;
for k = 1:N
  y = y + h*(y+2*x); % new value of y (computed from old x,y)
  x = x + h;         % new value of x
  plot(x, y, 'r*')
end

  % plot the exact solution (as a broken line) for comparison
x = linspace(x0, x0+N*h, 20);
y = exp(x) - 2*x -2;
plot(x,y)


  % B) input data h, xmax, x0, y0, f(x,y)
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
h = 0.1;           % step-size
xmax = 0.9;        % end-point of the interval
x0 = 0; y0 = -1;   % initial condition - starting point
plot(x0, y0, 'r*') % plot the starting point as a red asterisk
hold on            % continue on the same picture

  % explicit Euler
x = x0; y = y0;
while x < xmax
  y = y + h*(y+2*x);    % new value of y 
  x = x + h;            % new value of x
  plot(x, y, 'r*')
end


  % C) output data as vectors
  %%%%%%%%%%%%%%%%%%%%%%%%%%%
h = 0.1;             % step-size
xmax = 0.9;          % end-point of the interval
x(1) = 0; y(1) = -1; % initial condition - starting point

  % explicit Euler
k = 1;
while x(k) < xmax
  x(k+1) = x(k) + h;               % new value of x
  y(k+1) = y(k) + h*(y(k)+2*x(k)); % new value of y 
  k = k+1;
end

plot(x, y, 'r*')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% interpolation, approximation of given data-points
% by polynomials
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all
close all
N=9;

  % given data
x=linspace(0,4,N);
y=sin(x)-0.2*cos(4*x);

  % graph of the given data - black circles
plot(x,y,'ko','linewidth',2);
hold on

xx=linspace(0,4,40);  % mesh for graph of a polynomial

  % approximation by a polynomial using the least squares
  % try orders of the polynomial: 3, 6, 7, .. (up to N-1 ... interpol.) 
p = polyfit(x, y, 6); 
yy = polyval(p, xx);
plot(xx,yy,'m','linewidth',2);


  % ilustration of a typical behavior of polynomials of higher orders
y=[1 2 1 1 1 3 2 1 -1]; % for N=9; compare order 7 and 8
figure;   % new figure window
plot(x,y,'ko','linewidth',2);
hold on
p = polyfit(x, y, 7); 
yy = polyval(p, xx);
plot(xx,yy,'m','linewidth',2);  % magenta
p = polyfit(x, y, 8); 
yy = polyval(p, xx);
plot(xx,yy,'linewidth',2);      % blue

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%