Ok. I had some time to look at your surface more closely, It does not fit that PDE. It cannot do so. Now that I have your data...
S is a 45x84 array. Given your plot label,s, time moves from 0-83, and distance from 0-44. (Pick your own units for those variables.)
I''ve plotted S, as a function of distance at time 0. You want to use a diffusion model on this system, but have posed no forcing term in the equation of any form. Therefore, since this is a diffusion process, as time moves foward, this system will approach an AVERAGE level. Diffusion is an averaging thing. The level at a long time in the future, IF that model makes sense at all,
would be somethign near 1559.
But it CANNOT move to a much higher average level, unless your model includes a forcing term of some sort, where there is addtional "stuff" being injected over time. There is no such term in the model you posed.
So the line in red is the process after a relatively long time, the blue line is the process at time 0. Sorry, but this system is NOT a simple diffusion process. You CANNOT fit that model to this surface. The two are completely inconsistent.
Unless you have a better model for this process, trying to do the fit you asked to do would be meaningless here...
===================================
Edit:
I ran out of time last night to finish what I wanted to say. So now that we know the PDE that you posed CANNOT possibly solve the problem, is there a simple variation thereof that would do so? This now becomes entirely a problem in mathematical modeling, not for MATLAB, at least not immediately. Once we know what mathematics we need to use, then we can come back to MATLAB to write code.
Looking at the surface plot you show, we see the spatial variation is actually pretty small. Far more so than the temporal variation. So, temporarily, I'll assume that this is purely a temporal problem, so a 1-d ODE. What ODE has a solution that looks like an exponential rise to a constant asymptote? That one is easy.
dS/dt = -k*S + u
Here we don't know the rate constant. It might even be a function of x. But if we integrate that ODE, we would get a solution that behaves like your process.
dsolve(diff(S,t) == -k*S + u)
Here we would see a process that rises to an upper asymptote of u. The rate constant for the exponential is k. And C1 is an undetermined constant.
So now we can put this back into a PDE setting. It looks like after a long time, the process rises to an upper asymptote, one that is the same for all x. But the rate constant may vary. So I'll assume that k is an unknown function of x, though for a first pass, I might assume it is truly constant. Anyway, the PDE now becomes...
St = -k(x)*S + u + D*Sxx
Here St is the first partial of S with respect to t, and Sxx is the second partial with respect to x.
If k(x) is truly a constant, then the PDE becomes
St = -k*S + u + D*Sxx
This is still NOT something we can throw into PDEPE. But, since we have data, we can now use that data.
We can approximate those partial derivatives using finite differences. Then with k, u, and D as unknowns, the problem will become a simple linear least squares. In fact, we could then use backslash to estimate the unknown constants.
The next question is, what finite difference approximations would we employ? (By the way, all of this is enough that I could probably write a paper on the modeling, analysis and solution.) But your data is noisy. Enough so that a simple finite difference will be far more noisy. So I might use a trick from Savitsky-Golay now. For example, perhaps we might find a 7 point finite difference approximation to a series for the first derivative at any point, assuming the function is locally quadratic.
-0.10714 -0.071429 -0.035714 2.1371e-17 0.035714 0.071429 0.10714
The vector StFDA represents the coefficients of the best finite difference estimator of the first derivative of a vector, subject to noise, and assuming the process is locally quadratic. We can apply that finite difference approximation to the surface array using conv2.
Similarly, we can find another such approximation for the second partial derivative of the surface, though I'll need a factor of 2 to get the necessary coefficients.
0.11905 2.4547e-17 -0.071429 -0.095238 -0.071429 -3.2735e-18 0.11905
I know, you still don't see where this is going. But I do, and since this is an interesting problem in modeling, press forwards...
St_approx = conv2(S,StFDA','same');
Sxx_approx = conv2(S,SxxFDA,'same');
St_approx_trim = St_approx(3:end-2,3:end-2);
Sxx_approx_trim = Sxx_approx(3:end-2,3:end-2);
S_trim = S(3:end-2,3:end-2);
kuD = [-S_trim(:),ones(numel(S_trim),1),Sxx_approx_trim(:)]\-St_approx_trim(:)
So my best estimate of the rate constant k is 0.11611. The estimated upper asymptote is u/k,
and my estimate for D is 0.7523. Note that the upper asymptote is not too bad, compared to the data.
We see the upper asymptote does seem to pretty much pass through that final time point.
With a little more effort, I could have chosen a better form for k(x), if we wanted it to be non-constant. (k really does seem to vary with x, if I look at that surface.) But what I've already done is probably way more work than you ever thought was going to be necessary.
Finally, now that we know the estimates of k, u, and D, if we choose to provide intelligent boundary conditions and boundaty values on the solution, we could even throw this into PDEPE. The resulting prediction would now pretty nicely approximate the surface plot from your data.
