% discretization

function [yd,IS]=discry1(yd,fy,my,t,dis,m,Ify)

  yl=-1;                              % lower bound of output y
  yu=1;                              % upper bound of output y
  yr=yu-yl;
  y1=yl+dis*yr/m;
  y2=yl+dis*yr/m+yr/m;
  if fy(t)<=y1
    yd=1;
  elseif fy(t)<=y2
    yd=2;
  else 
    yd=3;
  end
  
  IS=(0.01/(my-1))*ones(my,my^2);
  IS(3,:)=0.99*ones(1,my^2);