function [PA,ID] = MyDisc(L,N)
% Adapted algorithm from Chebfun. Don't modify.
% L -- chebop, N = number of grid points
% PA -- corresponding matrix discritized, ID -- correponding identity on
% the same grid.
% in this code, L is chebop, A is corresponding linop.

prefs = cheboppref();

% Linearize, thereby obtaining linearity information, a LINOP, and an input of
% the correct dimensions to pass to N:
[A, ~, ~, ~] = linearize(L, L.init, [], true);

% Determine the discretization.
prefs = determineDiscretization(L, length(A.domain), prefs);

% Discretization type:
discType = prefs.discretization;

% Set up the discretization of A:
% Create a discretization object
discA = discType(A);

% Update the discretiztion dimension on unhappy pieces:
discA.dimension = repmat(N, 1, numel(discA.domain)-1);

if ( isempty(A.continuity) )
     % Apply continuity conditions:
     discA.source = deriveContinuity(discA.source);
end

% Discretize the operator (incl. constraints/continuity):
[PA, ~, ~, ~, PS] = matrix(discA);

% Identity operator
ID = [zeros(2,N+2); PS];
