clc;clear all;cellNum=10;
a = 10; %Width of slab (cm)
D =0.3733; %Diffusion coeffecient (cm)
sigA =0.0158; %Absorbption cross section (1/cm)
sigS =0.0158; %Scatt cross section (1/cm)
sigT =0.0158; %Total cross section (1/cm)
Source = 1; % source strength (neutrons/(cm3*s))
dx =a/cellNum; % Cell width
flux = sym('a', [1 cellNum]);
for i =1:cellNum
flux(i)=sym('a');
end
%Right Boundary Reflecting -[D d_flux/dx] = 0
B(1)=0;
A(1)=sigA;
%Left Boundary Condition
flux(cellNum)=0;
A(cellNum)=(2*D/dx^2+sigA);
%Middle Boundary conditions
for i=2:cellNum-1
C(i)=-D/dx^2;
B(i)=-D/dx^2;
A(i)=2*D/dx^2+sigA;
end
M = diag(A.*ones(1,cellNum)) + diag(B.*ones(1,cellNum-1),1) + diag(C.*ones(1,cellNum-1),-1);
Q = flux.*M;
Z=sum(Q,2);
for i=1:cellNum-1;
flux(cellNum-i)=solve(Z(cellNum-i)==1,flux(cellNum-i));
end
flux(1:end)
round(flux,4)
what do you mean you cant call the values in matrix? i think its working.
if you have an inital flux and dont want to change your inital flux, then you just need to change last few lines as
for i=1:cellNum-1;
fluxnew(cellNum-i)=solve(Z(cellNum-i)==1,flux(cellNum-i));
end
