function plot2ndODEs(fnct, indptvar, Iy, It, c1, c2, sp)
% Dr. M. Sakalli
% note dsolve expects D#y = ... Diff y and equality in text.. 
% in fact there are two solutions, I don't knw how 
% matlab determines to present solutions in array of 2, 
% or a superposed solution..

% Pr1_2ndOrderDE('D2y=-.0*Dy-(2*pi*2)^2*y+24*cos((2*pi*2)^2*t/2)', 't', [-.015 .015], [-.0 .5], 0.01, 0.01, 1)

% Pr2_2ndOrderDE('D2y=-.01*Dy-.02*y', 't', [-5 5], [1, 40], 2, 4, 5)
% fnct function, 
% indtvar;
% c1 c2 arbitrary variables. 
% It, Iy: Intervals. 
% choose sp an integer btw 1 =< sp<=5 ie..
% try someother functions, implicit solutions are okey.. 
% in the case of implicit solutions you need to give mesh x and y
% in the case of explicit solutions.. as below.. 

% sp for the sparcity of the arrow.. not that good..

sol=dsolve(fnct, indptvar),
soll=subs(sol,{'C1', 'C2'},{c1,c2}),
scalefactor=1;%=1/25*(t(50)-t(1))

if (~diff(It)|~diff(Iy))
	disp(sprintf('One of the intervals you entered is zero!, quiting'));
	return;
end

num=50;
xh=(It(2)-It(1))/num;      
yh=(Iy(2)-Iy(1))/num;
[x,y] = meshgrid(It(1):xh:It(2),Iy(1):yh:Iy(2));

%if solution is implicit this won't work..
% why need to take diff og solution.. 
slopes_=subs(diff(soll),'t',x);

rr=sqrt(xh.^2+slopes_.^2);
dxx=(1./rr)*scalefactor;
dyy=(slopes_./rr)*scalefactor;

% for the sparcity of the arrow.. not that good..
if ((sp-round(sp))|(sp<1)|(sp>5)) 
	sp=2;
end
quiver(x(1:sp:end, 1:sp:end),y(1:sp:end, 1:sp:end),dxx(1:sp:end, 1:sp:end),dyy(1:sp:end, 1:sp:end),.7, 'r');        % Here is the direction field
hold on
plot(It(1):xh:It(2),subs(soll,'t',It(1):xh:It(2)), '-*')
axis tight
hold off

return
