Nonlinear Regression on a Function Containing a Summation

I have been trying to fit the following dataset:

x=[0 0.0005 0.001 0.0015 0.002 0.003 0.004 0.005 0.0075 0.01 0.015 0.02 0.04 0.06 0.07 0.08 0.1 0.115];
y=[0 14229 28508 42737 56916 85374 99852 103846 114830 120821 133802 143787 168001 178486 179734 178985 169998 139793];

In the following equation, which contains a summation:

The following "myfun" function was defined as:

function output = myfun(coef,x)
%UNTITLED2 Summary of this function goes here
%   Detailed explanation goes here
nTerms=length(x); % The length of each function
yy=zeros(1,nTerms);
for i=1:1:length(x)
    P=zeros(1,nTerms);% Initiallize the sum to zero
    for k=1:1:nTerms
        n=k-1;
    P(k)=(((-1).^n)*(coef(1).^n)*(x(i).^(coef(4)*(n+1)-n*coef(3))))/((coef(2).^(n+1))*gamma(coef(4)*(n+1)-n*coef(3)+1));  % computing the sum
    end
    yy(i)=sum(P);
end
output=yy;
end

Initial value

x0=rand(1,4);

Solve nonlinear curve-fitting (data-fitting) problems in least-squares sense:

beta=lsqcurvefit(@(coef,x) myfun(coef,x),x0,x,y);

Finally, plot, calculate R^2 and RMSE

f=myfun(beta,x);
plot(x,y,'.')
plot(x,f,'-')

My problem is: Where do I have to apply the restrictions? Is the function well defined? How can I find better initial values, since each time the result changes?I tried other methods such as: fminsearch, lsqnonlin, etc.I hope you can help me identify my mistake.