how to simulate from a bivarite poisson(lamda1, lamda2, lamda3) using matlab

hi there please i need help in trying to simulate a poisson distribution with three lamda.i know for a single poisson we: clc clear all format compact rng(970,'v4') X =[]; U =[]; p =[0.3329 0.3662 0.2014 0.0738 0.0203 0.0045 0.0008 0.0001] for i = (1:10000) r =rand; if r>0 & r<=0.3329 x=0; X=[X,x]; U=[U,r]; else if r>0.3329 & r<=0.6991 x=1; X=[X,x]; U=[U,r]; else if r>0.6991 & r<=0.9005