There are many issues with your question. The most important one is this: your question is not urgent.
After that, you do an element-wise comparison for the input char arrays, instead of using strcmp (or switch-case). The root of your problem is that you try to write the offset elements from the diagonal, but don't check if those are actually inside the matrix.
You also seem to think Matlab is 0-indexed. The first array element in Matlab is not retrieved with a 0, but with a 1. I tried to fix the explicit mentions, but you should check them.
Additionally, you had a typo in your SORSC function: a lower case s instead of upper case, and a typo where you used len instead of length.
The code below still doesn't work as intended, but I can't correct that, because your loops are unclear what they should do. You overwrite the result every iteration, but it doesn't even start, because 1:x doesn't return a valid array to iterate over.
%####********Analytical Solutions *********#####
%#### This is simply plugging in the results from problem 1
%### Constants
sigmaa=0.005;
D=0.6;
L = (D/sigmaa)^0.5;
T=45;
So=10^8;
%### Slab Problem
x=linspace(0,45,10000);
U = (D/(2*L)-1/4)/(D/(2*L)+1/4);
C = (-So*T/sigmaa)*(T/4+D)/(-(D/(2*L)-1/4)^2/(D/(2*L)+1/4)*exp(-T/L)+(D/(2*L)+1/4)*exp(T/L));
A = U*C;
flux_slab=[];
for i = 1:x
ans1 = A*exp(-i/L)+C*exp(i/L)+(So*i^2)/sigmaa;
flux_slab = ans1;
end
%### Cylindrical Problem
x2=linspace(0,45/2,10000);
flux_cyl=[];
for i = 1:x2
ans1 = -So/(2*sigmaa*(0.5*besseli(0,T/(2*L))+(D/L)*besseli(1,T/(2*L))))*besseli(0,i/L)+(So/sigmaa);
flux_cyl = ans1;
end
%#### Sphere Problem
flux_sph=[];
for i = 1:x2
flux_sph = (-So*sinh(i/L)/(i*4*sigmaa*(sinh(T/(2*L))/(2*T)+(2*D/L)*(T/2*cosh(T/(2*L))-L*sinh(T/(2*L)))/(T^2))))+So/sigmaa;
end
flux_comp_slab_10 = Diffusion('slab', 45, 10, 0.005, 0.6, 10^8);
flux_comp_slab_100 = Diffusion('slab', 45, 100, 0.005, 0.6, 10^8);
flux_comp_slab_1000 = Diffusion('slab', 45, 1000, 0.005, 0.6, 10^8);
flux_comp_slab_10000 = Diffusion('slab', 45, 10000, 0.005, 0.6, 10^8);
flux_comp_cyl_10 = Diffusion('cylinder', 45/2, 10, 0.005, 0.6, 10^8);
flux_comp_cyl_100 = Diffusion('cylinder', 45/2, 100, 0.005, 0.6, 10^8);
flux_comp_cyl_1000 = Diffusion('cylinder', 45/2, 1000, 0.005, 0.6, 10^8);
flux_comp_cyl_10000 = Diffusion('cylinder', 45/2, 10000, 0.005, 0.6, 10^8);
flux_comp_sph_10 = Diffusion('sphere', 45/2, 10, 0.005, 0.6, 10^8);
flux_comp_sph_100 = Diffusion('sphere', 45/2, 100, 0.005, 0.6, 10^8);
flux_comp_sph_1000 = Diffusion('sphere', 45/2, 1000, 0.005, 0.6, 10^8);
flux_comp_sph_10000 = Diffusion('sphere', 45/2, 10000, 0.005, 0.6, 10^8);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PLOTTING OPTIONS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x3=linspace(0,45/2,10);
x4=linspace(0,45/2,100);
x5=linspace(0,45/2,1000);
x6=linspace(0,45/2,10000);
x7=linspace(0,45/2,100000);
figure(1);
plot(x, flux_slab);
hold on
plot(x3*2, flux_comp_slab_10);
plot(x4*2, flux_comp_slab_100);
plot(x5*2, flux_comp_slab_1000);
plot(x6*2, flux_comp_slab_10000);
hold off
figure(2);
plot(x, flux_cyl);
hold on
plot(x3*2, flux_comp_cyl_10);
plot(x4*2, flux_comp_cyl_100);
plot(x5*2, flux_comp_cyl_1000);
plot(x6*2, flux_comp_cyl_10000);
hold off
figure(2);
plot(x, flux_sph);
hold on
plot(x3*2, flux_comp_sph_10);
plot(x4*2, flux_comp_sph_100);
plot(x5*2, flux_comp_sph_1000);
plot(x6*2, flux_comp_sph_10000);
hold off
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% FUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DEFINITIONS
function A = SOR(x, So)
S = zeros(length(x),1); % BLANK MATRIX
for i = 1:(length(S)-1)
S(i)=So*(x(i)^2);
end
A=S;
return
end
function B = SORSC(x, So)
S = zeros(length(x),1); % BLANK MATRIX
for i = 1:(length(S)-1)
S(i)=So;
end
B=S;
return
end
function C = Diffusion(geometry, Size, N, SigmaA, D, So)
x = linspace(0,Size,N) ; % Set the size to iterate and number of nodes
delta = (x(2)-x(1)); % Distance between nodes
L = zeros(length(x), length(x)); % Start with a blank matrix
% GEOMETRY TYPE
switch geometry
case 'slab' % Depending on geometry make tridiagonal terms
S = SOR(x, So) / D; % Recieve source term from source def
for i = 1:(length(x)-1) % Iterate over tri diagonal
L(i,i) = 2/(delta^2)+SigmaA/D; % Diagonal Terms
if (i+1)<size(L,2)
L(i,i+1) = -1/(delta^2); % Above diagonal term
end
if (i-1)>0
L(i,i-1) = -1/(delta^2); % Below diagonal term
end
end
% Set boundary conditions at ends (Vacuum conditions)
L(1,1) = 0.5*(0.5-D/delta);
L(1,2) = D/(2*delta);
L(length(x),length(x)) = 0.5*(0.5-D/delta);
L(length(x), length(x)-1) = D/(2*delta);
case 'cylinder'
S = SORSC(x, So) / D ; % Recieve source term from source def
for i = 1:(length(x)-1) % Iterate over tri diagonal
L(i,i) = (x(i+1)/(x(i)*(delta^2))+1/(delta^2)+SigmaA/D); % Diagonal Terms
if (i+1)<size(L,2)
L(i,i+1) = -x(i+1)/(x(i)*delta^2); % Above diagonal term
end
if (i-1)>0
L(i,i-1) = -1/(delta^2); % Below diagonal term
end
end
% Set boundary conditions at ends (Vacuum conditions 2/ symmetry)
L(1,1) = -1/delta;
L(1,2) = 1/delta;
L(length(x),length(x)) = D/(2*delta);
L(length(x), length(x)-1) = 0.5*(0.5-D/delta);
case 'sphere'
S = SORSC(x, So) / D; % Recieve source term from source def
for i = 1:(length(x)-1) % Iterate over tri diagonal
% Diagonal Terms
L(i,i) = (x(i+1)^2)/((x(i)^2)*(delta^2)) + 1/(delta^2) + SigmaA/D;
if (i+1)<size(L,2)
L(i,i+1) = -x(i+1)^2 / (x(i)^2 * delta^2); % Above diagonal term
end
if (i-1)>0
L(i,i-1) = -1/(delta^2); % Below diagonal term
end
end
% Set boundary conditions at ends (Vacuum conditions w/ symmetry)
L(1,1) = -1/delta;
L(1,2) = 1/delta;
L(length(x),length(x)) = D/(2*delta);
L(length(x), length(x)-1) = 0.5*(0.5-D/delta);
otherwise
error('Not a valid geometry option');
end
flux = linsolve(L,S);
C=flux;
return
end
