function [R]=pkares(k,h,obsdata)
	k=abs(k);
	k1=[k 0]'; k2=[0 k]';
	f=@(c,h) k1./(k1+h).^2-k2./(k2+h).^2;

	c0=[1-k(1)/(k(1)+h(1)) k(1)/(k(1)+h(1)) zeros(1,rows(k1)-2)]; #yield initial guesses for concs
	C=lsode(f,c0,h);
	R=obsdata-C*pinv(C)*obsdata;
	R=sum(sum(sqrt(R.^2)));
endfunction
