function y=getpulse(x0,omega,T1,T2,mz,ka,t)

#Don't edit this part
	omegarf=omega(1);
	omega1=omega(2);
	omega0a=omega(3);
	omega0b=omega(4);

	T1a=T1(1);
	T1b=T1(2);
	
	T2a=T2(1);
	T2b=T2(2);
	
	mza0=mz(1);
	mzb0=mz(2);
	pab=mza0/mzb0;
	x=x0;

#You can edit this part if you know what you're doing
#To change rate law, edit e.g. -ka*x(1)+ka*pab*x(4) on each line

#Bloch equations follow here
mag= @(x,t) 	[ -(omegarf-omega0a)*x(2)-(x(1)/T2a)-ka*x(1)+ka*pab*x(4),
				(omegarf-omega0a)*x(1)+x(3)*omega1-(x(2)/T2a)-ka*x(2)+ka*pab*x(5),
				-x(2)*omega1-((x(3)-mza0)/T1a)-ka*x(3)+ka*pab*x(6),
				-(omegarf-omega0b)*x(5)-(x(4)/T2b)+ka*x(1)-ka*pab*x(4),
				(omegarf-omega0b)*x(4)+x(6)*omega1-(x(5)/T2b)+ka*x(2)-ka*pab*x(5),
				-x(5)*omega1-((x(6)-mzb0)/T1b)+ka*x(3)-ka*pab*x(6)];

	y=lsode(mag,x0,t);


endfunction
