efficient storing of matlab shape functions -

i'm trying find efficient possible way store , call matlab shape-functions. have interval x=linspace(0,20) , position-vector

count = 10; i=1:count     pos(i)=rand()*length(x); end 

and now, want put on every position pos(j) shape functions gauss-kernels compact support or hat-functions or similar (it should possible change prototype function easily). support of function controlled so-called smoothing length h. constructed function in .m-file (e.g. cubic spline)

function y = w_shape(x,pos,h)  l=length(x); y=zeros(1,l); if (h>0)     i=1:l         if (-h <= x(i)-pos && x(i)-pos < -h/2)             y(i) = (x(i)-pos+h)^2;         elseif (-h/2 <= x(i)-pos && x(i)-pos <= h/2)             y(i) = -(x(i)-pos)^2 + h^2/2;         elseif (h/2 < x(i)-pos && x(i)-pos <=h)             y(i) = (x(i)-pos-h)^2;         end     end else     error('h must positive') end 

and construct functions on interval x like

w = zeros(count,length(x)); j=1:count     w(j,:)=w_shape(x,pos(j),h); end 

so far good, when make x=linspace(0,20,10000) , count=1000, takes computer (intel core-i7) several minutes calculate whole stuff. since should kind of pde-solver, procedure has done in every time-step (under particular circumstances). think problem is, use x argument function-call , store every function instead of store 1 , shift it, matlab-knowledge not good, advices? fyi: need integral of areas, 2 or more function-supports intersect...and when i'm done in 1d, wanna 2d-functions, has efficient anyways

one initial vectorization remove loop in w_shape function:

for i=1:l     if (-h <= x(i)-pos && x(i)-pos < -h/2)         y(i) = (x(i)-pos+h)^2;     elseif (-h/2 <= x(i)-pos && x(i)-pos <= h/2)         y(i) = -(x(i)-pos)^2 + h^2/2;     elseif (h/2 < x(i)-pos && x(i)-pos <=h)         y(i) = (x(i)-pos-h)^2;     end end 

could become

xmpos=x-pos; % compute once , store instead of computing numerous times inds1=(-h <= xmpos) & (xmpos < -h/2); y(inds1) = (xmpos(inds1)+h).^2; inds2=(-h/2 < xmpos) & (xmpos <= h/2); y(inds2) = -(xmpos(inds2).^2 + h^2/2; inds3=(h/2 < xmpos) & (xmpos <=h); y(inds3) = (xmpos(inds3)-h).^2; 

there better optimisations this.

edit: forgot mention, should use profiler find out slow!

profile on run_code profile viewer