%Differentiation function
%the first parameter represents the time
%the second parameter represents the variable to be differentiated
%note that the two vectors must have the same size for the 
%function to work

%the function returns the derivative of y at the specified times t
%both inputs must be column vectors!

function ydot=derivative(t,y)
%You do not need to undersnatd how the function
%performs the differentiation.  If you are interested
%we recommend taking 16.901 next term
[s1,s2]=size(t);
[s3,s4]=size(y);
%check if the arguments are appopriate
if(s2~=1)
  error('time is not a column vector');
end;
if(s4~=1)
  error('y is not a column vector');
end;
if(s1~=s3)
  error('the time vector and the y vector differ in size');
end;
%t is a column vector containing the times
%y is a column vector containing the function values
%note: y is a vector, not a matrix!!!!
%Central Differencing Scheme
t1=circshift(t,+1);
t2=circshift(t,-1);
z1=circshift(y,+1);
z2=circshift(y,-1);
%derivative calculation
ydot=(z1-z2)./(t1-t2);
%make the lastone the same as the next to the last one
%because the derivative at the last point cannot be computed
[l,m]=size(ydot);
ydot(l,1)=ydot(l-1,1);
%same is true for the first point!
ydot(1,1)=ydot(2,1);