% This program will allow a user to pick various starting points (Initial Conditions)
% for solving of a system of differential equations.  The function will call
% g.m and use the canned differential equation solver ODE45 provided with Matlab.

clear all; close all;

N=5;
time=1;
xmin=-10;
xmax=10;
ymin=-10;
ymax=10;
zmin=-10;
zmax=50;
axis([xmin, xmax, ymin,ymax,zmin, zmax]);
%subplot(2,1,2), axis([xmin, xmax, ymin,ymax,zmin, zmax]);

hold on;

x(1)=-1;
y(1)=-1;
for i=2:N+1
   x(i)=x(i-1)+2/N;
   y(i)=y(i-1)+2/N;
end
for i=1:N+1
  for j=1:N+1
     zt(i,j)=1;
     end
  end
  % top and bottom
  
  mesh(x,y,zt);
  mesh(x,y,-zt);
for i=1:N+1
   xl(i)=1;
   for j=1:N+1
      zl(i,j)=y(j);
   end
end
 mesh(xl,y,zl);
 mesh(-xl,y,zl);
     
for i=1:N+1   
   for j=1:N+1
   [T,Y]=ode45('g',[0,time],[x(i),y(j),zt(i,j)]);
   plot3(Y(:,1),Y(:,2),Y(:,3),'r')
   nn=length(Y(:,1));
   plot3(Y(nn,1),Y(nn,2),Y(nn,3),'k*')
   [T,Y]=ode45('g',[0,time],[x(i),y(j),-zt(i,j)]);
    plot3(Y(:,1),Y(:,2),Y(:,3),'b')
   nn=length(Y(:,1));
    plot3(Y(nn,1),Y(nn,2),Y(nn,3),'k*')
   [T,Y]=ode45('g',[0,time],[xl(i),y(i),zl(i,j)]);
    plot3(Y(:,1),Y(:,2),Y(:,3),'g')
    nn=length(Y(:,1));
    plot3(Y(nn,1),Y(nn,2),Y(nn,3),'k*')
   [T,Y]=ode45('g',[0,time],[-xl(i),y(i),zl(i,j)]);
    plot3(Y(:,1),Y(:,2),Y(:,3),'g')
    nn=length(Y(:,1));
    plot3(Y(nn,1),Y(nn,2),Y(nn,3),'k*')

   
end

end