how to use lsqcurvefit for system of PDEs?is it possible?

I have a system of coupled PDEs (with a matrix as input)that ultimately gives me a matrix. I have written a code for solving this system numerically and after a specific time (user defined) it gives me the output matrix, this matrix should be matched with real data(a real image (matrix)). I have some coefficients that need to be optimized so that my output (predicted image) be matched as much as possible to the real image. I tried to use lsqcurvefit but it doesn't work. here is the code I've written: for n=2:nt

    for i=2:nx-1
        for j=2:ny-1
            u0(1,j)=0;
            u0(nx,j)=0;
            u0(i,1)=0;
            u0(i,ny)=0;

            %--------------------------------------------------------------
            u(i,j)=d1*(u0(i+1,j)+u0(i-1,j))+d2*(u0(i,j+1)+u0(i,j-1))+(1-2*(d1+d2))*u0(i,j)+(rho0(i,j)+((rhoT(i,j)-rho0(i,j))*t)/365)*dt*(u0(i,j)-(u0(i,j).^2));
            ICVFn(i,j)=u(i,j);
             if n==300
                 u(i,j)=inline('0.0206*(u0(i+1,j)+u0(i-1,j))+0.0206*(u0(i,j+1)+u0(i,j-1))+(1-2*(2*0.0206))*u0(i,j)+(rho0(i,j)+((rhoT(i,j)-rho0(i,j))*t)/365)*dt*(u0(i,j)-(u0(i,j).^2))','x','uo');
                 x=lsqcurvefit(u(i,j),x0,u0(i,j),ICVFT);
             end
        end
    end
    u;
    n;       
    t=t+dt;
    u0=ICVFn;

the error is: Error using inlineeval (line 15) Error in inline expression ==> 0.0206*(u0(i+1,j)+u0(i-1,j))+0.0206*(u0(i,j+1)+u0(i,j-1))+(1-2*(2*0.0206))*u0(i,j)+(rho0(i,j)+((rhoT(i,j)-rho0(i,j))*t)/365)*dt*(u0(i,j)-(u0(i,j).^2)) Undefined function 'u0' for input arguments of type 'double'.

Could you please help me to solve this? Does lsqcurvefit work for such a problem? if not, what can I use instead?