Nevermind. I figured it out.
Modelfun1 = @(beta1,xx)(beta1(1)-beta1(2).*(exp(xx./(beta1(4).*beta1(3)))-1));
Modelfun2 = @(beta2,xx)(beta2(1)-beta2(2).*(exp(xx./(beta2(4).*beta2(3)))-1));
beta1 = [I1(1) 1E-6 26 0.8];
beta2 = [I2(1) 1E-6 26 0.92];
Poly1 = nlinfit(U1,I1,Modelfun2,beta1);
Poly2 = nlinfit(U2,I2,Modelfun2,beta2);
Poly3 = polyfit(U1,P1,4); Poly4 = polyfit(U2,P2,4);
Trend1 = Modelfun1(Poly1,U); Trend2 = Modelfun2(Poly2,U);
Trend3 = polyval(Poly3,U); Trend4 = polyval(Poly4,U);
hold on
plot(U,Trend1,'k','LineWidth',1.5);
hold on
plot(U,Trend2,'k','LineWidth',1.5);
hold on
...
I had to design the function in Modelfun and estimate the coefficiants in beta.