matlab拟合空间曲线。散点拟合三维曲线

问题：已知三组数据，就是直角坐标系中的x，y，z的三个数组，要求在matlab中拟合出这条空间曲线。 帖子源数据：
a=

[22.77,17.13,14.37,12.98,12.34,12.12,12.11,12.18,12.27,12.32,12.32,12.27,12.18,12.11,12.12,12.34,12.98,14.37,17.13,22.77];
b=[23.75,21.25,18.75,16.25,13.75,11.25,8.75,6.25,3.75,1.25,-1.25,-3.75,-6.25,-8.75,-11.25,-13.25,-16.25,-18.75,-21.25,-23.75];
c=[0,3.34,6.59,9.44,11.82,13.74,15.21,16.29,16.99,17.33,17.33,16.99,16.29,15.21,13.74,11.82,9.44,6.59,3.34,0];

关于这个问题在论坛里被提及了好多次，貌似还没有最佳答案。在matlab中找了很多函数去实现，都不是很理想，也有人提到1st0pt曲线拟合软件，下面提供一种可行的解决方案：
step1：三组数据导入1stOpt（这是一个1分钟入手软件，即开即用，哈哈），程序-快速三维公式拟合，优化算法麦夸特法（Levenberg-Marquardt, LM）。选取相关系数为1的函数结构，同时还会生成多项式系数p。输出结果：
Function: z = p1+p2x+p3x^2+p4*y+p5*y2+p6*y^3
Algorithms: 麦夸特法(Levenberg-Marquardt) + 通用全局优化法
p1 21.414214338773
p2 -0.472199954302881
p3 0.0119553684329481
p4 0.00491434397562593
p5 -0.0298899195279345
p6 -9.69888338793035E-6

step2：在Matlab中系数p继续做优化处理， pp = lsqcurvefit(@fun,x0,xdata,ydata)。将上述p作为lsqcurvefit函数的x0预测值，fun为z方程式。求得优化后的系数pp。具体 lsqcurvefit见help文档

原始曲线和拟合曲线