How to fit non-linear matlab data?

Question

I'm working on a problem in scientific computing namely fitting data to this equation

$c(z) = 4800 + p_1 + p_2 \cdot z/1000 + p_3 \cdot e^{ -p4 \cdot z/1000}$

The data is in a background question that I asked previously to prepared for the problem

acoustics under water

Now to split my work into smaller units, I should first try and fit the data. In the answer to the question is says that I can use lsqnonlin but it doesn't say how. I look in the manual for lsqnonlin and it says that the syntax is

x = lsqnonlin(fun,x0)

but I don't fully understand how to use it. Ideally I would read my data from a file, can you tell me how to do it? The data is

z c(z)
0 5050
500 4980
1000 4930
1500 4890
2000 4870
2500 4865
3000 4860
3500 4860
4000 4865
5000 4875
6000 4885
7000 4905
8000 4920
9000 4935
10000 4950
11000 4970
12000 4990

and the problem specification is problem U8 in this text.

Update

Running the code from the answer gives the following result for me after 400 iterations:

ans = 3.6277 14.8277 22.2631 -0.1114

So this means p1=3.6277, p2=14.8277, p3=22.2631, p4=-0.1114

Is that correct?

Update 2

>> z=0:10:12000;
>> fittedfun=4800 + 3.6277+14.8277/1000*z+22.2631*exp(0.1114/1000*z);
>> plot(z, fittedfun);

It looks rather linear, so I'm not totally certain that I'm seeing what I expect.

Update 3

Comparing the fitted function with the actual data the approximation doesn't look very good so it leads me to think that there is something wrong with my method or the equation.

>> plot(z, fittedfun);
>> z=0:500:12000;
>> fittedfun=4800 + 3.6277+14.8277/1000*z+22.2631*exp(0.1114/1000*z);
>> plot(z, fittedfun);
>> data=[5050 4980 4930 4890 4870 4865 4860 4860 4865 4875 4885 4905 4920 4935 4950 4970 4990];
>> z=[0:500:4000 5000:1000:12000];
>> zu=z

zu =

  Columns 1 through 6

           0         500        1000        1500        2000        2500

  Columns 7 through 12

        3000        3500        4000        5000        6000        7000

  Columns 13 through 17

        8000        9000       10000       11000       12000

>> plot(data,zu)
>> plot(zu,data)
>> hold on;
>> plot(zu,data)
>> fittedfun=4800 + 3.6277+14.8277/1000*z+22.2631*exp(0.1114/1000*z);
>> z=[0:500:4000 5000:1000:12000];
>> plot(z, fittedfun);
>> 

Answer 1

First, collapse the Function to only necessary parameters (all constants are "eaten up"): $$ c(z) = p_1 + p_2 z + p_3 \exp(p_4 z) = p_1 + p_2 z + p_3 \tilde{z}^{p_4}$$ Now iterate over $p_4$ in plausible values and do a 2D-polyfit with $x=z, y=\exp(z)^{p_4}$, see http://www.mathworks.com/matlabcentral/newsreader/view_thread/17164 (google matlab 2d polyfit). Finding $p_4$ can be tricky, though.

For the nonlinear solver, use an error measure (i.e. fun = @(pVec)norm(data-c(p1, p2, p3, p4)) where pVec(i)=p_i).

EDIT:
The following MATLAB snippet did the work for me, though used maxFunEvals (400) so no idea how good the approximation is:
z=[0:500:4000 5000:1000:12000];
data=[5050 4980 4930 4890 4870 4865 4860 4860 4865 4875 4885 4905 4920 4935 4950 4970 4990];
fun=@(p)4800 + p(1)+p(2)/1000*z+p(3)*exp(p(4)/1000*z)-data;
x0=[0 1 1 -1]; % Guess
lsqnonlin(fun,x0) Another option would be to exploit $z=0$ to reduce the dimension if we force $c(0) = 5050$, but only if that is justified.

EDIT 2

The following script gives excellent results: z=[0:500:4000 5000:1000:12000]; data=[5050 4980 4930 4890 4870 4865 4860 4860 4865 4875 4885 4905 4920 4935 4950 4970 4990]; fun=@(p)(4800 + p(1))*ones(size(z)) +p(2)/1000*z+p(3)*exp(p(4)/1000*z)-data; x0=[0 0 0 -1]; % Guess opt = optimoptions('lsqnonlin', 'MaxFunEvals',1000); p=lsqnonlin(fun,x0,[],[],opt) fitf=@(t)(4800 + p(1))*ones(size(t)) + p(2)/1000*t+ p(3)*exp(p(4)/1000*t); tt=linspace(0,12000,1000); plot(z,data,'r-',tt,fitf(tt),'b-');

The only change was another guess in $x0$.