How can I adjust the output of my pdepe function for nlinfit?

Hi!

I have the following system of coupled PDEs which I have solved using the pdepe solver. I have compared it to an analytical solution and verified it. Now I am trying to fit some experimental data to this system. I am using the nlinfit function to do achieve this.

After that I calculate the flux using pdeval to use in below expression. I am using the final row of the solution BTW so basically .

And thus finally I get my

Now I guess the problem here is that I am unable to compute the necessary number of solutions for or I_num as it is referred to in the code.

As I am getting the following error.

Error using nlinfit (line 219)
MODELFUN must be a function that returns a vector of fitted values the same size as Y (6-by-1).  The model function you provided returned a result that was 1-by-1.

One common reason for a size mismatch is using matrix operators (*, /, ^) in your function instead of the corresponding elementwise operators (.*, ./, .^).

Error in PSw_EK_v4_Caller (line 26)
[beeta] = nlinfit(alphaa,I_exp,@PSt_EK_v4,beta0); 

I am calling the pdepe using the following script.

I_exp      = [2.5e-03; 
              2.8e-04; 
              2.5e-03; 
              2.8e-04; 
              2.5e-03; 
              2.8e-04];

alphaa     = [5.00000000000000e-06,3.46000000000000e-06,3.54444000000000e-07;
              5.00000000000000e-06,3.46000000000000e-06,3.54444000000000e-07;
              5.00000000000000e-06,3.46000000000000e-06,3.54444000000000e-07;
              5.00000000000000e-06,3.46000000000000e-06,3.54444000000000e-07;
              5.00000000000000e-06,3.46000000000000e-06,3.54444000000000e-07;
              5.00000000000000e-06,3.46000000000000e-06,3.54444000000000e-07];

% wrapper    = @(beeta, alphaa)arrayfun(@(si)PSt_EK_v4(beeta,si), alphaa);
beta0      = [5e-06; 0; 2.81E-06;  800;    1e-05; 1e05];             % starting point 
[beeta] = nlinfit(alphaa,I_exp,@PSt_EK_v4,beta0);