% A script to run continuous time systems
%
% The system is of the form x' = f(x)
%
% The function "f" is defined in a file called "rhs.m"
%
global FNAME;
global GNAME;
disp('Gradient system');




if (length(FNAME)==0), FNAME='frhs'; end;
entry = input(['Enter name of the function (', FNAME ,')-> '],'s');
if (length(entry) == 0), entry = FNAME; end;
FNAME = entry;



if (length(GNAME)==0), GNAME='grhs'; end;
entry = input(['Enter name of the gradient (', GNAME ,')-> '],'s');
if (length(entry) == 0), entry = GNAME; end;
GNAME = entry;




x0 = input('Please enter the initial vector -> ');
[r,c] = size(x0);
if (r==1), x0 = x0'; end;
t0 = input('Please enter the starting time (default 0) -> ');
if (isempty(t0)), t0=0; end;
t1 = input('Please enter the ending time (default 1) -> ');
if (isempty(t1)), t1=1; end;

disp('Working.');

[t,x] = ode45('gpatch', t0,t1,x0);


disp('Computing graph.');


xmin = min(x(:,1));
xmax = max(x(:,1));

ymin = min(x(:,2));
ymax = max(x(:,2));

steps = 20;

xlist = [xmin : (xmax-xmin)/steps : xmax ];
ylist = [ymin : (ymax-ymin)/steps : ymax ];

[xx,yy] = meshdom(xlist,ylist);

zz = zeros(steps+1,steps+1);
dx = zz;
dy = zz;

for p = 1:steps+1, for q = 1:steps+1
	zz(p,q) = feval(FNAME,[xx(p,q);yy(p,q)]);
	g = feval(GNAME,[xx(p,q);yy(p,q)]);
	dx(p,q) = g(1);
	dy(p,q) = g(2);
end;end;

zmin = min(min(zz));
zmin = min(zmin,0);
zmax = max(max(zz));
zmax = max(zmax,0);


n = length(t);
z = zeros(n,1);

for k=1:n;
	z(k) = feval(FNAME,x(k,:)');
end;

disp('Plotting.')


mesh(xx,yy,zz);
hold on;
extra = 0;
plot3(x(:,1),x(:,2),z + extra);
plot(x(:,1),x(:,2));  % the "shadow"
quiver(xx,yy,-dx,-dy);
contour(xx,yy,zz,'-g');
axis([xmin,xmax,ymin,ymax,zmin,zmax]);
grid;

hold off;