function syspde_para
clear all;
close all;
% n_eq = number of equations
n_eq = 4; %INPUT

% m_int = number of spatial intervals
m_int = 20; %INPUT
xo = 0.0;
xf = 1.0;
n_int = 1000;
to = 0.0;
tf = 200.0;
%  BC is a vector of boundary condition type definitions
%  0 = Dirichlet
%  1 = Neumann
BC = [0 1; 0 1; 0 1; 0 1;]; %INPUT

%  code must create files
% function y_out = syspde_boundary_conditions(k_eq,k_bc,x,t,y1,i_x, dx);
% function y_out = syspde_initial_conditions(k_eq,x);





%
%  create spatial grid
%
%  mx = number of spatial gridpoints
%  (this grid has imaginary points for Neumann BCs)
%  (it still works if the problem has only Dirichlet BCs)
mx = m_int + 3;
dx = (xf-xo)/m_int;
dxi = 1/dx;
xgrid = [(xo-dx):dx:(xf+dx)]';
%  create temporal grid
%  nt = number of temporal gridpoints
nt = n_int + 1;
dt = (tf-to)/n_int;
tgrid = [(to):dt:(tf)]';
%
%  dimension solution matrix
%
y = zeros(n_eq,nt,mx);
%
%  define some limitations on indices for where solution is not given by ICs and
%  where the derivatives need to be calculated using finite difference techniques
%
for k_eq = 1:1:n_eq
   for i_BC = 1:1:2
      if (i_BC == 1)
   		if (BC(k_eq, i_BC) == 0) 
      		i_x_v(k_eq,i_BC) = 3;
   		else
      		i_x_v(k_eq,i_BC) = 2;
         end
      else
   		if (BC(k_eq, i_BC) == 0) 
      		i_x_v(k_eq,i_BC) = mx-2;
   		else
      		i_x_v(k_eq,i_BC) = mx-1;
         end
      end
   end
end
%i_x_v
%i_x_d
%
%  input initial conditions into solution matrix
%
j_t = 1;
for k_eq = 1:1:n_eq
   for i_x = 2:1:mx-1
      x = xgrid(i_x);
      y(k_eq,j_t,i_x) = syspde_initial_conditions(k_eq,x);
   end
end
%
%  input Dirichlet Boundary conditions into solution matrix for all times
%
y1 = zeros(n_eq,mx);
for k_eq = 1:1:n_eq
   for k_bc = 1:1:2
   	if (BC(k_eq,k_bc) == 0) 
         if (k_bc == 1) 
            i_x = 2;
         else
            i_x = mx - 1;
         end
			x = xgrid(i_x);
         for j_t = 1:1:nt
            t = tgrid(j_t);
            for m_eq = 1:1:n_eq
               y1(m_eq,:) = y(m_eq,j_t,:);
            end
%            fprintf(1,'D: k_eq= %i kbc= %i i_x = %i x = %f t = %f  \n',k_eq, k_bc,i_x,x,t)
            y(k_eq,j_t,i_x) = syspde_boundary_conditions(k_eq,k_bc,x,t,y1,i_x,dx);
         end
      end
   end
end
%
%  input initial Neumann Boundary conditions into solution matrix
%
j_t = 1;
t = tgrid(j_t);
for k_eq = 1:1:n_eq
	y1(k_eq,:) = y(k_eq,j_t,:);
end
for k_eq = 1:1:n_eq
   for k_bc = 1:1:2
   	if (BC(k_eq,k_bc) == 1) 
         if (k_bc == 1) 
            i_x = 1;
            i_x2 = 2;
            i_x3 = 3;
         else
            i_x = mx;
            i_x2 = mx-1;
            i_x3 = mx-2;
         end
         x = xgrid(i_x);
%         fprintf(1,'N: k_eq= %i kbc= %i i_x = %i x = %f t = %f  \n',k_eq, k_bc,i_x,x,t)
			y(k_eq,j_t,i_x) = syspde_boundary_conditions(k_eq,k_bc,x,t,y1,i_x,dx);
      end
   end
end
%
%  call Classical second order Runge-Kutta method 
%
dydt_1 = zeros(n_eq,mx);
dydt_2 = zeros(n_eq,mx);
for j_t = 2:1:nt
   %  obtain first half of approximation to derivative 
   t = tgrid(j_t-1);
   for k_eq = 1:1:n_eq
		y1(k_eq,:) = y(k_eq,j_t-1,:);
	end
   dydt_1 = syspde_para_input(n_eq,i_x_v,y1,t,dxi,mx);
   %  obtain second half of approximation to derivative 
   t = t + dt;
      %  first obtain intermediate values of y in interior nodes
   for k_eq = 1:1:n_eq
      for i_x = i_x_v(k_eq,1):1:i_x_v(k_eq,2)
         y1(k_eq,i_x) = y1(k_eq,i_x) + dt*dydt_1(k_eq,i_x);
      end
   end
   %  second obtain intermediate values of y from Dirichlet BCs
   for k_eq = 1:1:n_eq
   	for k_bc = 1:1:2
   		if (BC(k_eq,k_bc) == 0) 
          	if (k_bc == 1) 
            	i_x = 2;
         	else
            	i_x = mx-1;
	         end
         	x = xgrid(i_x);
            y1(k_eq,i_x) = y(k_eq,j_t,i_x);
         end
      end
   end
   % third obtain intermediate values of y from Neumann BC's
	t = tgrid(j_t);
	for k_eq = 1:1:n_eq
   	for k_bc = 1:1:2
   		if (BC(k_eq,k_bc) == 1) 
         	if (k_bc == 1) 
            	i_x = 1;
         	else
            	i_x = mx;
	         end
         	x = xgrid(i_x);
				y1(k_eq,i_x) = syspde_boundary_conditions(k_eq,k_bc,x,t,y1,i_x,dx);
      	end
   	end
	end
   dydt_2 = syspde_para_input(n_eq,i_x_v,y1,t,dxi,mx);
   for k_eq = 1:1:n_eq
      for i_x = 1:1:mx
         y(k_eq,j_t,i_x) = y(k_eq,j_t-1,i_x) + 0.5*dt*(dydt_1(k_eq,i_x) + dydt_2(k_eq,i_x) );
      end
   end
   % compute new values of unknown from Neumann BC's
	t = tgrid(j_t);
	for k_eq = 1:1:n_eq
   	for k_bc = 1:1:2
   		if (BC(k_eq,k_bc) == 1) 
         	if (k_bc == 1) 
            	i_x = 1;
         	else
            	i_x = mx;
	         end
         	x = xgrid(i_x);
				y(k_eq,j_t,i_x) = syspde_boundary_conditions(k_eq,k_bc,x,t,y1,i_x,dx);
      	end
   	end
	end
end
%
%  write solution
%
ix = 12
it = 501
tgrid(it)
xgrid(ix)
y(:,it,ix)
%
%  plot solution
%
% define plot/movie parameters
T=y;
%
% normalize plots by inlet values
%
Cscale = T(1,1,2);
  T(1,:,:) = T(1,:,:)/Cscale;
  T(2,:,:) = T(2,:,:)/Cscale;
  T(3,:,:) = T(3,:,:)/T(3,1,2);
  T(4,:,:) = T(4,:,:)/T(4,1,2);
%
nframe = 1000;
if (n_int > nframe)
	nskip = round(n_int/nframe);
else
   nskip = 1;
end
lplot = 1;
fps = 3;
%  define plot line colors (for up to 10 eqns)
plc{1} = 'k-';
plc{2} = 'r-';
plc{3} = 'b-';
plc{4} = 'g-';
plc{5} = 'm-';
plc{6} = 'k:';
plc{7} = 'r:';
plc{8} = 'b:';
plc{9} = 'g:';
plc{10} = 'm:';
if (lplot == 1)
   fpsinv = 1.0/fps;
	newplot;
	Tmax=max(max(max(T)));
   Tmin=min(min(min(T)));
   xtext = xo + 0.1*(xf-xo);
   ytext = Tmin +0.9*(Tmax-Tmin);
   for j_t=1:nskip:nt
      for k_eq = 1:1:n_eq
         for i_x = 2:1:mx-1
            vector_plot(i_x) = T(k_eq,j_t,i_x);
         end
         plot(xgrid(2:mx-1) , vector_plot(2:mx-1),plc{k_eq} );
         hold on
      end
      axis([xo xf Tmin Tmax]);
      XLABEL('position along pipe (m)') 
      YLABEL('Temperature (K)') 
      temps = char(strcat('time = ', num2str(tgrid(j_t)), ' sec '));
      text(xtext, ytext,temps)
      hold off
      pause(fpsinv);
	end
end

function y_out = syspde_boundary_conditions(k_eq,k_bc,x,t,y1,i_x, dx);
if (k_bc == 1)
	if (k_eq == 1) 
   	y_out = 10;
	elseif (k_eq == 2)
		y_out = 0;
	elseif (k_eq == 3)
		y_out = 45.6;
	else
		y_out = 300.0;
%      y_out = y1(k_eq, i_x+2)- 2*dx*(-1.0);
   end
else
	if (k_eq == 1) 
		y_out = y1(k_eq, i_x-2)- 2*dx*(0.0);
	elseif (k_eq == 2)
		y_out = y1(k_eq, i_x-2)- 2*dx*(0.0);
	elseif (k_eq == 3)
		y_out = y1(k_eq, i_x-2)- 2*dx*(0.0);
	else
		y_out = y1(k_eq, i_x-2)- 2*dx*(0.0);
   end
end
%fprintf(1,'BC: k_eq= %i kbc= %i i_x = %i x = %f t = %f yout = %f \n',k_eq, k_bc,i_x,x,t,y_out)

function y_out = syspde_initial_conditions(k_eq,x);
if (k_eq == 1) 
   y_out = 0;
elseif (k_eq == 2)
	y_out = 0;
elseif (k_eq == 3)
	y_out = 55.6;
else
   y_out =  300;
end