%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
%  Solve the differential equation 
%
%  y' = F(t,y)
%
%  with the initial condition y(Ti) = Yi
%
%  on the time interval [Ti,Tf] with N timesteps, using an unspecified
%  mystery method.
% 
%
%
%  Input:
%     F - the right-hand side function of (t,y)
%    Ti - an intial time
%    Tf - a final time
%    Yi - the initial y-value
%     N - number of timesteps
%
%  Returns:
%
%    t - an array of N+1 equally spaced times, starting at Ti and ending at Tf
%    y - an array of approximate solution values at each of the t's
%
%  Example:
%
%  F = @(t,y) y 
%  [t,y] = MysterySolver(F,0,1,1,100)
%  plot(t,y) 

function [t,y]=MysterySolver(f,Ti,Tf,Yi,N)
  t = zeros(1,N+1);
  y = zeros(1,N+1);  
  h = (Tf-Ti)/N;

  t(1) = Ti;
  y(1) = Yi;
  for( k=1:N )
    tk=t(k); yk=y(k);
    k1 = h*f(tk,yk);
    k2 = h*f(tk+h/2,yk+k1/2);
    k3 = h*f(tk+h/2,yk+k2/2);
    k4 = h*f(tk+h,yk+k3);
    t(k+1) = Ti+k*h;
    y(k+1) = y(k) + (k1/6+k2/3+k3/3+k4/6);
    
  end
end