How to solve Diffusivity equation with sensitivity analysis using single code for changing variable values

Hi all,

I am solving this 1D-diffusivity equation problem with following code. The code works perfectly fine. But I am looking to solve it for multiple values of various parameters. I am trying to create multiple solutions using same code rather than changing values in the code every time.

e.g. Starting with Ct value (initially 1e-5 in the code), I want to solve the equation for values of 1e-5 and 1e-6. I can plug the values separately and get answers. But I want to know how to solve it with for both values in one code.

I tried to solve for Ct = [1e-5, 1e-6]. But it got messed up at Amat(i, i) = 2+alpha;, throwing error that Subscripted assignment dimension mismatch (Line 38).

%% 1D Pressure Diffusivity Equation
%Grid specification
nx = 10;
dx = 10;            % ft
%% Time specifications
dt = 1;         %day
ti = 0;
tf = 5;
%% Boundary conditions
P1 = 1000;          % psia
P10 = 500;           % psia
%% Initial conditions
Pi = 1000;          % psia
%% Formation properties
phi = 0.2;
mu = 100;           % cp
Ct = 1e-5;          % 1/psia
k = 10;            % md
%% Array allocations
Amat = zeros(10, 10);
bvector = ones(10,1);
P_nplus1 = ones(10,1);      %Pressure at future condition
P_n = ones(nx,1);           %Pressure at current condition
%% Computation
alpha = phi*mu*Ct/(0.00633*k)*(dx)^2/dt;
%% Assign BC and IC
P_n(1,1) = P1;
P_n(nx,1) = P10;
P_n(2:(nx-1),1) = Pi;
P_save = ones(nx, (tf-ti/dt));
%% Array elements assignment
Amat(1,1) = 1;
Amat(10,10) = 1;
for i=2:9
    Amat(i, i) = 2+alpha;
    Amat(i, i-1) = -1;
    Amat(i, i+1) = -1;
    bvector(i,1) = alpha*P_n(i,1);
end
bvector(1,1)= P1;
bvector(nx,1) = P10;
%% Time proceeding
t = ti;
j = 0;
while t < tf
    j = j+1;
    for i = 2:(nx-1)
        bvector(i,1) = alpha*P_n(i,1);
    end
    P_nplus1 = Amat\bvector;
        P_save(:,j) = P_nplus1;
        P_n = P_nplus1;
        t = t + dt;

end
%% Matrix Solution
P_nplus1 = Amat\bvector;
%% Plots for sensitivity analysis
plot(P_save)