%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
%  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,x,y]=RK4(f,Ti,Tf,xi,yi,h)
  N = ceil((Tf-Ti)/h);
  t = zeros(1,N+1);
  x = zeros(1,N+1);
  y = zeros(1,N+1);  
  h = (Tf-Ti)/N;

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