%exercise 7.11
%define constants
global k


vf=50;
vs=5;
t=0;
u=[2000 0]';
U=[];
T=[];

%free fall
h=0.1;
%k=(m*g)/vf^2;
k=(75*9.81)/vf^2;
number_step=30/h;
for step=1:number_step
   f1=f711(t,u);
   f2=f711(t+h/2,u+h*f1/2);
   f3=f711(t+h/2,u+h*f2/2);
   f4=f711(t+h,u+h*f3);
   u=u+(h/6)*(f1+2*f2+2*f3+f4);
   U=[U; u'];
   T=[T; t];
   t=t+h;
end
disp('h(3)=');
disp(u(1));

%k=(m*g)/vs^2;
k=(75*9.81)/vs^2;

h=0.001;
number_step=2/h;
for step=1:number_step
   f1=f711(t,u);
   f2=f711(t+h/2,u+h*f1/2);
   f3=f711(t+h/2,u+h*f2/2);
   f4=f711(t+h,u+h*f3);
   u=u+(h/6)*(f1+2*f2+2*f3+f4);
   U=[U; u'];
   T=[T; t];
   t=t+h;
end
 
 h=0.1;
while u(1)> 0 
   f1=f711(t,u);
   f2=f711(t+h/2,u+h*f1/2);
   f3=f711(t+h/2,u+h*f2/2);
   f4=f711(t+h,u+h*f3);
   u=u+(h/6)*(f1+2*f2+2*f3+f4);
   
   U=[U; u'];
   T=[T ;t];
   t=t+h;
   
end
% % % %
disp('h(t) på marken');
disp(-u(2));
disp('Tid total');
disp(t);
Tend=length(T);
subplot(2,1,1), plot(T ,U(:,1), [0 180] ,[0 0], '--'), title('hÃ¶jden');
subplot(2,1,2), plot(T ,-U(:,2), [0 180] ,[0 0], '--'), title('hastighet');
