%% Stable
f = @(t,y) -y-sin(t)+cos(t);
tv = 0:.1:10; plot(tv,cos(tv),'b'); hold on % plot exact solution
[ts,ys] = Euler(f,[0,10],1,50); % use 50 steps of size 10/50
plot(ts,ys,'r.-');

%% Unstable
f = @(t,y) y-sin(t)-cos(t);
tv = 0:.1:5; plot(tv,cos(tv),'b'); hold on % plot exact solution
[ts,ys] = Euler(f,[0,5],1,50); % use 50 steps of size 5/50
plot(ts,ys,'r.-');

%% Stiff ODE

f = @(t,y) y.^2 - y.^3;
df = @(t,y) 2*y - 3*y.^2;%% Logistic

delta = .001;
h = 4;
N = ceil(2/(delta*h));
[ts,ys] = Euler(f,[0,2/delta],delta,N);
[t1s,y1s] = BEuler(f,df,[0,2/delta],delta,N);
plot(ts,ys,'o-r'); hold on;
plot(t1s,y1s,'o-b');
hold off;

axis([0 2/delta 0 1.2]);

%% Stiff ODE with MATLAB solvers

f = @(t,y) y.^2 - y.^3;
delta = .001;
[ts,ys] = ode45(f,[0,2/delta],.01);
[t1s,y1s] = ode15s(f,[0,2/delta],.01);
[t2s,y2s] = ode23(f,[0,2/delta],.01);
%plot(ts,ys,'o-r'); hold on;
%plot(t1s,y1s,'o-b');
plot(t2s,y2s,'o-g');
hold off;

axis([0 2/delta 0 1.2]);


%% Functions

function [ts,ys] = Euler(f,I,y0,n)
   ys = zeros(0,n);
   ys(1)=y0;
   T = I(2);
   h = T/n;
   ts = linspace(0,T,n);
   for i=1:n-1
       ys(i+1) = ys(i) + h*f(ts(i),ys(i));
   end
end

function x = Newton(f,df,x,max_iter,max_size)
      sqrteps = sqrt(eps(2));
      steps = 0;
      err = 3*sqrteps;
      while err > sqrteps
         h = f(x)./df(x);
         err = abs(h);
         x = x - h;
         steps = steps +1;
         if (steps > max_iter) || (abs(x) > max_size)
             err= 0;
         end
      end
end

function [ts,ys] = BEuler(f,df,I,y0,n)
   ys = zeros(0,n);
   ys(1)=y0;
   T = I(2);
   h = T/n;
   ts = linspace(0,T,n);
   g = @(t,y1,y0) y1 - h*f(t,y1) - y0;
   dg = @(t,y1) 1 - h*df(t,y1);
   for i=1:n-1
       gi = @(y1) g(ts(i),y1,ys(i));
       dgi = @(y1) dg(ts(i),y1);
       ys(i+1) = Newton(gi,dgi,ys(i),20,10);
   end
end