Vectorized Relaxation Algorithm For Solving Linear Partial Differential Equations

Utilizes a vectorized relaxation algorithm to solve linear partial differential equations.

This program utilizes a vectorized relaxation algorithm in order to solve linear partial differential equations. Please submit any suggestions you have for improving the program (I would eventually like to include a convergence constant in order to implement the SOR method).
Example:
%This is an example of how to use the program for a parallel plate capacitor with random point charges in the grid space:
%sets a plot constant to plot data train every 100 loop counts
plotk=100;
%Here are the grid dimensions (note that the grid need not be square)
L_xspan=150;
L_yspan=150;

meshdim=1; %set the unit height of the grid points
tol=10^(-7); %tolerance

%places the plates in the center of the grid with an x length 5
plate1=L_xspan/2-5;
plate2=L_yspan/2+5;

%separates the plates by 10 grid units
position_plate1=L_yspan/2-5:L_yspan/2+5;
position_plate2=L_yspan/2-5:L_yspan/2+5;

%creates the walls of the grid and capacitor plates
Bound_0_x={1:L_xspan;1:L_xspan;1;L_xspan;plate1;plate2};
Bound_0_y={1;L_yspan;1:L_yspan;1:L_yspan;position_plate1;position_plate2};

%sets the voltage of either plate
v1=175;
v2=-175;

%sets the values of the grid walls at zero and inserts the capacitor plate voltages
Bound_0_val=[0;0;0;0;v1;v2];

%generates 5 random gridpoints
Random_gridpoints=round(100*rand(2,5));

%generates 5 random potentials
Random_potentials=round(100*rand(1,5));

%inserts initial values for the random potentials that were generated
init_x={Random_gridpoints(1,1);Random_gridpoints(1,2);Random_gridpoints(1,3);Random_gridpoints(1,4);Random_gridpoints(1,5)};
init_y={Random_gridpoints(2,1);Random_gridpoints(2,2);Random_gridpoints(2,3);Random_gridpoints(2,4);Random_gridpoints(2,5)};
init_val=[Random_potentials(1,1);-Random_potentials(1,2);Random_potentials(1,3);Random_potentials(1,4);-Random_potentials(1,5)];

%calls the function
[array_f,X,Y] = FDRelax(meshdim,L_xspan,L_yspan,Bound_0_x,...
Bound_0_y,Bound_0_val,init_x,init_y,...
init_val,tol,plotk);