Matlab - Asistente tecnico

Hola, estoy teniendo un problemita con la convergencia del Metodo de Newthon Rhapson para sistemas de ecuaciones lineales. Lo que trato de obtener es apartir de una ecuacion dendo no conosco 3 parametros, pero tengo muchos puntos como para comprobarlo. Los parametros a averiguar son k, a y b. D, E son constantes y dp, P y hp son variables de las cuales tengo varios puntos.
dp/D=k*(P*E*D^-2)^a * (hp/D)^b

A continuacion le copio el codigo:

clear,clc,cla,close all
format long
tol=0.01;
imax=200;
x=[0.1,0.1,0.1];
[xr,kp]=newtonsi(x,tol,imax)

-----------------------------------------------

function [xr,kp]=newtonsi(x,tol,imax)
format long
kp=1;
epi=1;
x1=x;
while norm(epi)>tol
x=x1;
fxn=f(x)
axn=jac(x)
epi=axnfxn';
x1=x-epi';
kp=kp+1
if kp>imax
disp('no converge')
break
end
end
xr=x1;

----------------------------------------------------

function y=f(x)
% funcion para utilizar con newtonsi.m
format long
E=200000;
p=429.005;
D=1.5875;
hp=0.04256;
dp=0.53642;

y(1)=x(1)*((p*E*D^(-2))^x(2))*(hp/D)^x(3)-dp/D

E=200000;
p=1030;
D=1.5875;
hp=0.0877;
dp=0.7596;

y(2)=x(1)*((p*E*D^(-2))^x(2))*(hp/D)^x(3)-dp/D

E=200000;
p=1900;
D=1.5875;
hp=0.18;
dp=1.10988;

y(3)=x(1)*((p*E*D^(-2))^x(2))*(hp/D)^x(3)-dp/D

----------------------------------------------------------

function df=jac(x)
% matriz jacobiana para usar con newtonsi.m
format long

E=200000;
p=429.005;
D=1.5875;
hp=0.04256;
dp=0.53642;

df(1,1)=(p*E/D^2)^x(2)*(hp/D)^x(3);
df(1,2)=x(1)*(p*E/D^2)^x(2)*log(p*E/D^2)*(hp/D)^x(3);
df(1,3)=x(1)*(p*E/D^2)^x(2)*(hp/D)^x(3)*log(hp/D);

E=200000;
p=1030;
D=1.5875;
hp=0.0877;
dp=0.7596;

df(2,1)=(p*E/D^2)^x(2)*(hp/D)^x(3);
df(2,2)=x(1)*(p*E/D^2)^x(2)*log(p*E/D^2)*(hp/D)^x(3);
df(2,3)=x(1)*(p*E/D^2)^x(2)*(hp/D)^x(3)*log(hp/D);

E=200000;
p=1900;
D=1.5875;
hp=0.18;
dp=1.10988;

df(3,1)=(p*E/D^2)^x(2)*(hp/D)^x(3);
df(3,2)=x(1)*(p*E/D^2)^x(2)*log(p*E/D^2)*(hp/D)^x(3);
df(3,3)=x(1)*(p*E/D^2)^x(2)*(hp/D)^x(3)*log(hp/D);