% Euler's method

clear

% y' = x^2 example
h=1;N=round(10/h); % step size h and number N of steps 
x=(0:N)'*h; % x_k are multiples of h
y=zeros(N+1,1);
yI=zeros(N+1,1);
y(1)=1; % initial value
yI(1)=1;

for k=1:N

    y(k+1)=y(k)+h*x(k)^2; % Euler's method

    yI(k+1)=yI(k)+h/2*(x(k)^2+x(k+1)^2); % improved Euler's
end

truesoln=x.^3/3+1;
figure;plot(x,truesoln,'k',x,y,'r',x,yI,'b')
legend('True','Euler','Improved')

 % y' = y example
h=1;N=round(10/h);
x=(0:N)'*h;
y=zeros(N+1,1);
yI=zeros(N+1,1);
y(1)=1; % initial value
yI(1)=1;
for k=1:N

    y(k+1)=y(k)+h*y(k); % Euler's method

    p=yI(k)+h*y(k); % predictor
    yI(k+1)=yI(k)+h/2*(y(k)+p); % corrector
end

truesoln=exp(x);
figure;plot(x,truesoln,'k',x,y,'r',x,yI,'b')
legend('True','Euler','Improved')