% Matlab-skripti: suuntak.m
% Kts. myös exasuuntak.m
% HA 18.9.2000
%

clf
h=1;a=-2;b=2;c=-2;d=2;   % Muutettavat parametrit
t=a:h:b;x=c:h:d;                          
% f=inline('t.*x','t','x')   % Muutettava funktio
%  Laajempi esimerkkidata: suuntakdata.m
% m=1;b=1; g=9.81;               % Teht. 1 LV 10.9.01
% f=inline('g-(b*x.^2)/m','t','x')   % inline-määrittely ei hyväksy globaaleja.

% m=1;b=1; g=9.81;               
  f=inline('9.81-(x.^2)','t','x')   %  Teht. 1 LV 10.9.01    
  
 
[T,X]=meshgrid(t,x);
Z=f(T,X)
flipud(Z)   % havainnollistukseen, tosi käytössä kommentti tai puolipiste
plot(T(:),X(:),'x')
hold on
U=ones(size(Z));
normit=sqrt(1+Z.^2);
quiver(T,X,0.5*h*U./normit,0.5*h*Z./normit,0)
axis equal
axis([a-1,b+1,c-1,d+1])
grid
title('x''=tx')
shg
%
% Jos haluat pelkän suuntakentän, lopeta tähän.
%
% Seuraavassa piirretään Matlabin numeerisen ODE-ratkaisijan laskemia
% ratkaisukäyräapproksimaatioita hiirellä valittavaan alkupisteeseen.
%
display('Valitse hiirellä alkupisteitä, lopuksi ENTER grafiikkaikkunassa')
[t,x]=ginput;
for i=1:length(t)
  [T,Y]=ode45(f,[t(i),b+0.5],x(i));
  plot(T,Y,'r')
  [T,Y]=ode45(f,[t(i),a-0.5],x(i));
  plot(T,Y,'r')
end;  
shg






