clc; clear; close all; % known parameter c=0.95; siga=0.02; sigt=siga/(1-c); D=1/(3*sigt); nu_sigf=0.0165; N=100; H=100; h=H/N; % calculate parameter sigd=D/h^2; ave_siga=siga; ave_nu_sigf=nu_sigf; alpha_l=10^7; %black alpha_r=alpha_l; %black % set matrix for M and F (linear system) M=zeros(N+1); F=zeros(N+1); M(1,1)=sigd+(alpha_l/h)+(0.5*siga); M(N+1,N+1)=M(1,1); F(1,1)=0.5*nu_sigf; F(N+1,N+1)=F(1,1); for i=2:N M(i-1,i)=-sigd; M(i,i-1)=-sigd; M(i,i)=2*sigd+ave_siga; F(i,i)=ave_nu_sigf; end %%% solve matrix %%% using eig and inv A=inv(M)*F; A [flux,k]=eig(A) flux_eig1=flux(:,1); k_eig1=k(1,1); flux_eig2=flux(:,2); k_eig2=k(2,2); flux_eig3=flux(:,3); k_eig3=k(3,3); k_eig=[k_eig1 k_eig2 k_eig3]