Fitting Numerical Solution of PDE to Given Data with Multiple Parameters

I'm trying to fit my solution to be as close as possible to given data for my PDE system. I've been using the following as a reference (https://www.mathworks.com/matlabcentral/answers/406938-fitting-the-numerical-solution-of-a-pde-with-one-parameter-to-some-data) which has been super helpful. I was able to reproduce the results from the example (shown below), and I've adapted the code to my PDE system (which I'm solving for 10 parameters instead of one).

But as I run my code, the parameters are not adjusting at all, and it produce the following result.

As the message said, I adjusted the step tolerance to be lower (10^(-10)) so that it may continue find the optimal parameters. But as its being solved, the parameters don't adjust from the initialized ones. Is there a possible explanation and a remedy for this issue?

(Note: I'll provide my code as a comment to this question post. But of course, I won't be able to provide the data I'm using to compare to train the optimization algorithim )

global exper 
exper=[-0.4,-112,-121,-128,-131,-131.1,-133.2,-133.1,-134.4,-134.7,-135.12,-135.7,-135.8,-135.1,-135.4,-135.8,-135.9,-134,-135.7,-135.1,-135.3,-135.4,-136,-139,-136]; 
% p=[0.218,0.52,0.0115,2.03,31.6]; 
% qr f a n ks 
p0=[0.218,0.52,0.0115,2.03,31.6]; 
lb=[0.01,0.4,0,1,15]; %lower bound 
ub=[Inf,Inf,10,10,100]; %upper bound 
options = optimoptions('lsqnonlin','Algorithm','trust-region-reflective','Display','iter'); 
[p,resnorm,residual,exitflag,output] = lsqnonlin(@richards1,p0,lb,ub,options) 
uu=richards1(p); 
global t 
plot(t,uu+exper') 
hold on 
plot(t,exper,'ko') 
function uu=richards1(p) 
% Solution of the Richards equation     
%    using MATLAB pdepe                   
% 
%   $Ekkehard Holzbecher  $Date: 2006/07/13 $ 
% 
% Soil data for Guelph Loam (Hornberger and Wiberg, 2005) 
% m-file based partially on a first version by Janek Greskowiak 
%-------------------------------------------------------------------------- 
L = 200;                  % length [L] 
s1 = 0.5;                 % infiltration velocity [L/T] 
s2 = 0;                   % bottom suction head [L] 
T = 4;                    % maximum time [T] 
% qr = 0.218;               % residual water content 
% f = 0.52;                 % porosity 
% a = 0.0115;               % van Genuchten parameter [1/L] 
% n = 2.03;                 % van Genuchten parameter 
% ks = 31.6;                % saturated conductivity [L/T] 
x = linspace(0,L,100); 
global t 
t = linspace(0,T,25); 
options=odeset('RelTol',1e-4,'AbsTol',1e-4,'NormControl','off','InitialStep',1e-7) 
u = pdepe(0,@unsatpde,@unsatic,@unsatbc,x,t,options,s1,s2,p); 
global exper 
uu=u(:,1,1)-exper'; 
% figure; 
% title('Richards Equation Numerical Solution, computed with 100 mesh points'); 
% subplot (1,3,1); 
% plot (x,u(1:length(t),:)); 
% xlabel('Depth [L]'); 
% ylabel('Pressure Head [L]'); 
% subplot (1,3,2); 
% plot (x,u(1:length(t),:)-(x'*ones(1,length(t)))'); 
% xlabel('Depth [L]'); 
% ylabel('Hydraulic Head [L]'); 
% for j=1:length(t) 
%     for i=1:length(x) 
%         [q(j,i),k(j,i),c(j,i)]=sedprop(u(j,i),p); 
%     end 
% end 
%   
% subplot (1,3,3); 
% plot (x,q(1:length(t),:)*100) 
% xlabel('Depth [L]'); 
% ylabel('Water Content [%]'); 
% ------------------------------------------------------------------------- 
function [c,f,s] = unsatpde(x,t,u,DuDx,s1,s2,p) 
[q,k,c] = sedprop(u,p); 
f = k.*DuDx-k; 
s = 0; 
end 
% ------------------------------------------------------------------------- 
function u0 = unsatic(x,s1,s2,p) 
u0 = -200+x; 
if x < 10 u0 = -0.5; end 
end 
% ------------------------------------------------------------------------- 
function [pl,ql,pr,qr] = unsatbc(xl,ul,xr,ur,t,s1,s2,p) 
pl = s1; 
ql = 1; 
pr = ur(1)-s2; 
qr = 0; 
end 
%------------------- soil hydraulic properties ---------------------------- 
function [q,k,c] = sedprop(u,p) 
qr = p(1);               % residual water content 
f = p(2);                 % porosity 
a = p(3);               % van Genuchten parameter [1/L] 
n = p(4);                 % van Genuchten parameter 
ks = p(5);                % saturated conductivity [L/T]   
m = 1-1/n; 
if u >= 0 
    c=1e-20; 
    k=ks; 
    q=f; 
else 
    q=qr+(f-qr)*(1+(-a*u)^n)^-m; 
    c=((f-qr)*n*m*a*(-a*u)^(n-1))/((1+(-a*u)^n)^(m+1))+1.e-20; 
    k=ks*((q-qr)/(f-qr))^0.5*(1-(1-((q-qr)/(f-qr))^(1/m))^m)^2; 
end 
end 
end 

[Thomas Rodriguez]

Initial Parameters are given as which are the same values called upon from (PDE_Coefficients(9)).

Code:

clc; clear; close all;

% Defining Actual Data:clc; clear; close all;

% Defining Actual Data:
global actual
Actual_1D = load('Actual_1D.mat', 'Actual_1D');
Actual_1D = Actual_1D.Actual_1D;
actual = Actual_1D;

% Defining Initial Parameters:
[D1,D2,B_e,B_s,B_a,B_v,theta,lambda,gamma,Omega,rho_e,rho_s,rho_a] = PDE_Coefficients(9);
% Storing Initial Parameters that are being varied:
p0 = [theta, lambda, gamma, B_e, B_v, D1, D2, rho_e, rho_s, rho_a];
% lower bound
lb=[10^(-4), 10^(-4), 10^(-4), 10^(-15), 10^(-20), 10^(-8), 10^(-8), 1000, 10^(6), 10^(4)];
% upper bound
ub=[1.0, 1.0, 1.0, 10^(-4), 10^(-4), 10^(-1), 10^(-1), 10^(4), 9*10^(8), 10^(6)];

options = optimoptions('lsqnonlin','Algorithm','trust-region-reflective','Display','iter'); % ,'StepTolerance',10^(-9),'MaxIterations',30 ,'FunctionTolerance',10^(-8)
[p,resnorm,residual,exitflag,output] = lsqnonlin(@PDE_Model,p0,lb,ub,options)

% Plotting:
global t x
global loc
loc = 12; % Specific Location being examined
UU = PDE_Model(p);
figure; hold on; grid on; box on;
plot(t,UU + actual(:,loc)) 
hold on;
plot(t,actual(:,loc),'ro')
xlabel('t')
ylabel('U')
title(['Population Density at ' num2str(x(loc)) '$$^\circ$$'], 'Interpreter','latex')
legend('Model','Actual Data','Location','best')


figure; hold on;
global Approx
imagesc(x,t,Approx); % Contour Plot
grid on; box on;
xlabel('Latitude$^\circ$',"Interpreter","latex","FontSize",12)
ylabel('Time (Weeks)',"Interpreter","latex","FontSize",12)
title('Approximated: $$\theta \lambda \hat{E}$$',"Interpreter","latex","FontSize",12)
axis tight
colormap jet
colorbar
hold off;

figure; hold on;
imagesc(x,t,actual); % Contour Plot
grid on; box on;
xlabel('Latitude$^\circ$',"Interpreter","latex","FontSize",12)
ylabel('Time (Weeks)',"Interpreter","latex","FontSize",12)
title('Actual: $$\theta \lambda E$$',"Interpreter","latex",'FontSize',12)
axis tight
colormap jet
colorbar
hold off;

% Printing Table:
% p = [theta, lambda, gamma, B_e, B_v, D1, D2, rho_e, rho_s, rho_a];
fprintf('\n')
fprintf('Optimal Parameter Values       \n')
fprintf('-----------------------\n')
fprintf('D1 %3.4e \n', p(6))
fprintf('D2 %3.4e \n', p(7))
fprintf('theta %3.4e \n', p(1))
fprintf('lambda %3.4e \n', p(2))
fprintf('gamma %3.4e \n', p(3))
fprintf('B_e %3.4e \n', p(4))
fprintf('B_v %3.4e \n', p(5))
fprintf('rho_e %3.4e \n', p(8))
fprintf('rho_s %3.4e \n', p(9))
fprintf('rho_a %3.4e \n', p(10))
fprintf('-----------------------\n')


function UU = PDE_Model(p)
% Defining Space and Time:
global A B
[A,B] = Loading_Data(); % Loading Dataframes
global x t t0 tf
% Space (x):
x = sort(A(:,1));
t0 = 9; tf = 37;
% Time (t):
t = [t0:1:tf]';
% Solving Equation:
fprintf('Solving PDE System: \n')
options = odeset('RelTol',1e-4,'AbsTol',1e-4,'NormControl','off','InitialStep',1e-7);
m = 0;
sol = pdepe(m,@pdefun,@pdeic,@pdebc,x,t,options,p);
% Solutions for the 6 Populations:
E = sol(:,:,2);
global actual loc
fprintf('Integrating Newly Infected Symptomatic Data, theta*lambda*E: \n')
N = length(t); theta = p(1); lambda = p(2);
global Approx
Approx = RK4(E,t0,tf,N,theta,lambda);
fprintf('Computing Error Matrix, UU: \n\n')
UU = Approx - actual;
% disp(UU)
UU = UU(:,loc);

% pdepe form:
% c(x,t,u,u_x)u_t = x^(-m)*(x^(m)*f(x,t,u,u_x))_x + s(x,t,u,u_x)
function [c,f,s] = pdefun(x,t,u,dudx,p)
% Parameters:
theta = p(1);
lambda = p(2);
gamma = p(3);
B_e = p(4);
B_v = p(5);
D1 = p(6);
D2 = p(7);
rho_e = p(8);
rho_s = p(9);
rho_a = p(10);
[B_s, B_a] = transmission_constant(t);
Omega = Omega_Surface_Trans(7*t);

% pdepe form:
% c(x,t,u,u_x)*u_t = x^(-m)*(x^(m)*f(x,t,u,u_x))_x + s(x,t,u,u_x)

% c(x,t,u,u_x):
% 6 populations for 6 system of eq'ns
c = ones([6 1]);

% f(x,t,u,u_x,u_y):
f = [D1; D1; D1; D1; D1; D2].*dudx;

% s(x,t,u,u_x):
F1 = -(B_e*u(2) + B_s*u(3) + B_a*u(4) + B_v*u(6)).*u(1);
F2 = (B_e*u(2) + B_s*u(3) + B_a*u(4) + B_v*u(6)).*u(1) - lambda*u(2);
F3 = theta*lambda*u(2) - gamma*u(3);
F4 = (1 - theta)*lambda*u(2) - gamma*u(4);
F5 = gamma*(u(3) + u(4));
F6 = rho_e*u(2) + rho_s*u(3) + rho_a*u(4) - Omega*u(6);
F = [F1;F2;F3;F4;F5;F6];
s = F;
end

% Initial Conditions:
function [u0] = pdeic(x,p)
% Input: scalar of space array (x) i.e. 1x1
% Initial Population:
[S0,E0,I_S0,I_A0,R0,V0] = Population(x);
% Storing Initial Pop. into u0:
u0 = [S0; E0; I_S0; I_A0; R0; V0];
end

% Establishing Boundary Conditions:
function [pl,ql,pr,qr] = pdebc(xl,ul,xr,ur,t,p)
% Standard Boundary Conditions Solver:
% p(x,t,u) + q(x,t)*f(x,t,u,u_x) = 0
% No Flux Boundary Conditions:
% Domain: (ymin = 32.55, ymax = 33.03)

% Boundary Conditions on the Left:
% S_x(32.55,t) = E_x(32.55,t) = I_s_x(32.55,t) = 
% I_a_x(32.55,t) = R_x(32.55,t) = V_x(32.55,t) = 0, etc.
pl = zeros([6 1]);
ql = ones([6 1]);

% Boundary Conditions on the Right:
% S_x(33.03,t) = E_x(33.03,t) = I_s_x(33.03,t) = 
% I_a_x(33.03,t) = R_x(33.03,t) = V_x(33.03,t) = 0, etc.
pr = pl;
qr = ql;
end

% Determing B_s and B_a:
function [B_s, B_a] = transmission_constant(t)
% B_s and B_a are defined as such in
% Angleica Bloomquist, Naveen Vaidya, et. al. ”Modeling Risk of Environmental Reservior of SARS-CoV-2”
% Input: t - time (scalar)
% t0 = 9; % initial time
% tf = 38; % final time
tol = 10^(-2); % tolerance condition
n = 1000; % number of points in stage arrays
% First Stage: 03/07/2020 - May 10, 2020 (~ 64 days ~ 9 weeks)
t1 = linspace(t0,18,n);
% Second Stage: May 11, 2020 - July 14,2020 (~ 64 days ~ 9 weeks)
t2 = linspace(18.01,28,n);
% Third Stage: July 15, 2020 - September 16, 2020 (~ 63 days ~ 9 weeks)
t3 = linspace(28.01,tf,n);
for i = 1:n
    if (t1(i) == t) | (abs(t - t1(i)) <= tol)
        % First Stage: 03/07/2020 - May 10, 2020 (~ 64 days ~ 9 weeks)
        B_s = 8.92*10^(-8);
        B_a = 2.67*10^(-8);
        break
    elseif (t2(i) == t) | (abs(t - t2(i)) <= tol)
        % Second Stage: May 11, 2020 - July 14,2020 (~ 64 days ~ 9 weeks)
        B_s = 3.32*10^(-8);
        B_a = 2.49*10^(-8);
        break
    elseif (t3(i) == t) | (abs(t - t3(i)) <= tol)
        % Third Stage: July 15, 2020 - September 16, 2020 (~ 63 days ~ 9 weeks)
        B_s = 1.64*10^(-8);
        B_a = 1.54*10^(-8);
        break
    else
        continue
    end
end
end

% Omega function:
function [omega] = Omega_Surface_Trans(t)
Final_tmp = Final_Temperature(t);
% Defining a and b values:
a = 0.094519;
b = -6.4829;
omega = a*Final_tmp + b;
% Note: a and b values are the coefficient averages over all surfaces
end

% Runge-Kutta (4th Order):
function Approx = RK4(E,t0,tf,N,theta,lambda)
h = (tf - t0)/N;
tn = [t0:h:tf]';
% Calculating Cumulative Data:
% f(t) = \theta*\lambda*E(t,x)
% E(t,x) assumed to be 29x30 size
% theta = 0.6;
% lambda = 0.21;
f = @(t,E) theta*lambda*E;
% Initial Condition: y(1) = y0
Y(1,:) = theta*lambda*E(1,:);
% Runge-Kutta Method:
for i = 1:N-1
    % k1 = f( t_n, y_n )
    k1 = f(tn(i), E(i,:));
    % k2 = f( t_n + 0.5*h, y_n + h*0.5*k1 )
    k2 = f(tn(i) + 0.5*h, E(i,:) + h*0.5*k1);
    % k3 = f( t_n + 0.5*h, y_n + h*0.5*k2 )
    k3 = f(tn(i) + 0.5*h, E(i,:) + h*0.5*k2);
    % k4 = f( t_n + h, y_n + h*k3 )
    k4 = f(tn(i) + 0.5*h, E(i,:) + h*k3);
    % f(t_n+1, y_n+1)
    Y(i+1,:) = Y(i,:) + (1/6)*h*(k1 + 2*k2 + 2*k3 + k4);
end
Approx = Y;
end

function [S0,E0,I_S0,I_A0,R0,V0] = Population(x)
% fprintf('x input: %2.4f \n', x)
[S,E,I_S,I_A,R,V] = Initial_Pop();
tol = 10^(-8); % Defined Tolerance:
% [A,~] = Loading_Data();
Lat = sort(A(:,1)); % Latitude Coordinates:
for i = 1:length(Lat)
    if (abs(Lat(i) - x) <= tol)
        S0 = S(i);
        E0 = E(i);
        I_S0 = I_S(i);
        I_A0 = I_A(i);
        R0 = R(i);
        V0 = V(i);
        break
    else
        continue
    end
end
end
end

% Calculating Initial Population Function:
function [S0,E0,I_S0,I_A0,R0,V0] = Initial_Pop()
global A B
Lat = A(:,1); Zip = string(A(:,4));
Pop = A(:,3); Case_Data = B(2,:); % Week 1: 03/1/2020 - 03/07/2020

% % Calculating Initial Population Distribution:
% Initial Susceptible Population Density:
[S0] = Initial_Population_Density(Lat,Pop,Zip,"Initial Susceptible $$S(0,X)$$ Population Density",0);
n = length(S0);
% Initial Exposed Population Density:
E0 = zeros(1,n);
% Initial Symptomatic Infected Population Density:
I_S0 = Initial_Population_Density(Lat,Case_Data,Zip,"Initial Susceptible $$S(0,X)$$ Population Density",0);
% Initial Asymptomatic Infected Population Density:
I_A0 = zeros(1,n);
% Initial Recovered Population Density:
R0 = zeros(1,n);
% Initial Environment Population Density:
V0 = zeros(1,n);
end

Final Temperature Function:

% Temperature Functions:
function [F] = Final_Temperature(x)
% Note: x = [0:(1/24):365]' gives the Net Temperature fluctation for One Year
f_s = Seasonal_Temp(x);
f_d = Dirunal_Temp(x);
F = f_s + f_d;
end

% Seasonal Temperature Fluctuation Function:
function [f_s] = Seasonal_Temp(x)
t_0 = 64.712;
epsilon_m = -6.9178;
tau_m = 365;
phi_m = 0.99605;
f_s = t_0 + epsilon_m*sin((2*pi/tau_m)*x + phi_m);
end

% Dirunal Temperature Fluctuation Function:
function [f_d] = Dirunal_Temp(x)
epsilon_d = -4.3372;
tau_d = 1;
phi_d = -5.0872;
m = 0;
f_d = m + epsilon_d*sin((2*pi/tau_d)*x + phi_d);
end