%Local and global errors for Euler's Method for Solving Initial Value Problems
%Use with ydot.m to evaluate rhs of differential equation
% Input: interval [a,b], initial value y0, step size h
% Output: time steps t, solution y
% Example usage: y=global_local2([0 1],1,0.1);
function [t,y]=global_local2(int,y0,h)
t(1)=int(1); y(1)=y0;
n=round((int(2)-int(1))/h);
tol=odeset('RelTol',1e-8);
[T Y] = ode45(@ydot,int,y0,tol);  
plot(T,Y,'Linewidth',2)          %plot "exact" solution Y
hold on
for i=1:n-1
  t(i+1)=t(i)+h;
  y(i+1)=eulerstep(t(i),y(i),h);
  for j=1:size(Y,1)
      if(T(j)<=t(i+1)) & (T(j+1)>t(i+1))  %find ts s.t. Y(ts)=y(i+1)
          ii=j;
      end
  end
  s=(t(i+1)-T(ii))/(T(ii+1)-T(ii));       
  ys=Y(ii)*(1-s)+Y(ii+1)*s;               % find Y(ts)
  [T1 Y1] = ode45(@ydot,[t(i+1) int(2)],y(i+1),tol);           
  plot(t,y,'r.')                                               % plot Euler values
  pause;
  line([t(i+1) t(i+1)], [y(i+1) ys],'Linewidth',2,'color','r') %draw local error line
  pause;
  plot(T1,Y1)            % plot next exact path starting in y(i+1)
  pause;
  T=T1;
  Y=Y1;
end
%hold off

function y=eulerstep(t,y,h)
%one step of Euler's Method
%Input: current time t, current value y,  stepsize h
%Output: approximate solution value at time t+h
y=y+h*ydot(t,y);

function z=ydot(t,y)
%z = t*y + t^3;
z=y*y*y -5*t;
