Add these lines after your existing code:

[x1Grid,x2Grid] = meshgrid(linspace(min(X(:,1)), max(X(:,1)), 100), ...
                          linspace(min(X(:,2)), max(X(:,2)), 100));
[~,scores] = predict(SVMModel, [x1Grid(:), x2Grid(:)]);
scores = reshape(scores, [size(x1Grid,1), size(x1Grid,2), 3]);

contour(x1Grid, x2Grid, scores(:,:,1) - max(scores(:,:,2:3),[],3), [0 0], 'k-');
contour(x1Grid, x2Grid, scores(:,:,2) - max(scores(:,:,[1 3]),[],3), [0 0], 'k--');
contour(x1Grid, x2Grid, scores(:,:,3) - max(scores(:,:,1:2),[],3), [0 0], 'k:');

The outcome is like: